{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "44394ee22333d032",
   "metadata": {},
   "source": "# Binary Companion "
  },
  {
   "cell_type": "code",
   "id": "357d57ec38373d8c",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-11-28T04:44:33.199378Z",
     "start_time": "2024-11-28T04:44:32.979020Z"
    }
   },
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams[\"font.size\"] = 20"
   ],
   "outputs": [],
   "execution_count": 1
  },
  {
   "cell_type": "code",
   "id": "692b3ee4af4f2d05",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-11-28T04:44:33.256279Z",
     "start_time": "2024-11-28T04:44:33.204268Z"
    }
   },
   "source": [
    "xx = np.linspace(0, 2*np.pi)\n",
    "yy = np.sin(xx)\n",
    "fig, ax = plt.subplots(1,1, figsize=(8, 5))\n",
    "ax.hlines(0, 0-0.5, 2*np.pi+0.5, linestyles=\"--\", color=\"k\")\n",
    "ax.plot(xx, yy, lw=3)\n",
    "\n",
    "ax.errorbar(.5*np.pi, 0.5, 0.5, 0, capsize=5, color=\"tab:orange\", )\n",
    "ax.text(1.7, 0.4, \"K [km/s]\", fontdict={\"color\":\"tab:orange\"})\n",
    "\n",
    "ax.errorbar(np.pi, 1.2, 0, np.pi, capsize=5, color=\"tab:orange\", )\n",
    "ax.text(2.8, 0.9, \"P [day]\", fontdict={\"color\":\"tab:orange\"})\n",
    "\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "ax.set_ylabel(\"Radial Velocity [km/s]\")\n",
    "ax.set_xlabel(\"Time\")\n",
    "ax.set_xlim(0-0.5, 2*np.pi+0.5)\n",
    "ax.set_ylim(-1.5, 1.5)\n"
   ],
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1.5, 1.5)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 800x500 with 1 Axes>"
      ],
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqEAAAG2CAYAAABPmV/gAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABu50lEQVR4nO3dd3gUVdsG8Ht2N70SkhBKGiWU0HsHEZGqCKF3RMWuiPraFf1UBBVUxALSO0gv0nsvSQgttCRACmkkIX135/tjScKkbtjNTrJ7/64rVzZnzsw8CbA8OWfOcwRRFEUQEREREZmQQu4AiIiIiMjyMAklIiIiIpNjEkpEREREJscklIiIiIhMjkkoEREREZkck1AiIiIiMjkmoURERERkciq5AygPrVaL6OhoODk5QRAEucMhIiIiokJEUURaWhpq1aoFhaLk8c4qlYRGR0fD29tb7jCIiIiIqAx37txBnTp1SjxepZJQJycnALpvytnZWeZoiIiIiKiw1NRUeHt75+dtJalSSWjeFLyzszOTUCIiIqJKrKxHJ7kwiYiIiIhMjkkoEREREZkck1AiIiIiMjkmoURERERkckxCiYiIiMjkmIQSERERkckxCSUiIiIik2MSSkREREQmxySUiIiIiEyOSSgRERERmRyTUCIiIiIyOSahRERERGRyTEKJiIiIyOSYhBIRERGRyTEJJSIiIiKTYxJKRERERCbHJJSIiIiITI5JKBERERGZHJNQIiIiIjI5JqFEREREZHJMQomIiIjI5JiEEhEREZHJMQklIiIiIpNjEkpEREREJscklIiIiIhMjkkoEREREZkck1AiIiIiMjkmoURERERkckxCiYiIiMjkmIRWhLRY4MB3us9EREREhjLD3IJJaEVIiwUOfW9Wf1GIiIhIRmaYWzAJJSIiIiKTYxJKRERERCbHJJSIiIiITE4ldwBmTZ0J5KTLHQURERFVdepMuSMwOiahFemfvnJHQERERFQpcTqeiIiIiEyOI6EVafIuwKu53FEQUWUQFwbEXyv/eR4NgRpNjR8PEVUtsaFmN8PKJLQiqewAawe5oyCiysC7g+6DiOhJqOzkjsDoOB1PRERERCbHJJSIiIiITI5JKBERERGZHJNQIiIiIjI5vRYmTZ48uaLjgCAIWLhwYYXfxyScvIAe/9N9JiIiIjKUGeYWgiiKYlmdFAoFBEGosCBEUYQgCNBoNKX2S01NhYuLC1JSUuDs7Fxh8RARERHRk9E3XytXiaaaNWvCysrK4OAel5ubi+joaKNek4iIiIgqt3Ilobt370aTJk2MGkBYWBiaN2dBdyIiIiJLIvvCpIqc5iciIiKiykn2JJSIiIiILI9e0/G3b98GANSuXdvoATRs2DD/+kRERERkGfRKQn19fSsuAJWqQq9PRERERJUPp+OJiIiIyOSYhBIRERGRyZWrRNOT2Lp1K9auXYuEhAT4+/tjypQpaN26dUXfloio8vvSpWibb1dg0vYnv+bPzYCUKKDFaOCF+U9+nYq0aAAQebRo+5cppo+FiGRjUBJ64MABjBgxAra2tggNDYWrq6vk+GeffYZvv/1W0rZgwQIsXLgQ48aNM+TWREQV5/YRYMnA4o+p7AAHd8CrORD4gu5DWeG/zxMRmR2D3jl37NiBhIQEvPDCC0US0NDQUHz77bfI2xW0WrVqSE5OhlqtxiuvvIJu3brBz8/PkNsTEZmeOhNIuaP7uLYdOPk7MGo14FTjya/Z9kWg3RTda2t748RZmQ2eB+Rk6F6fWQCcXShvPEQkC4OeCT169CgEQUDv3r2LHJs/fz5EUUS1atVw7tw5JCYm4vTp03Bzc0N2djb++OMPQ25NRGQabV8EXj1R8PHiXqDfLMDVR3c8+jywehTw6BfuJ+LgAdRoovuo5meUsCu1an4F36+Dh9zREJFMDEpCY2JiAACBgYFFjm3btg2CIOCNN95Aq1atAABt27bFG2+8AVEUsXfvXkNuTURkGo8niDWaAN7tgA4vA68cBtzq6vrcOweE75I3TiKiKsagJDQ+Ph4AikzF37x5E/fu3QMAvPDCC5Jj3bp1y+9DRFRl2VUDuk4r+PoGf7EmIioPg54JzXveMyVFuqLxyJEjAAAXFxe0bNlScqx69eoAgIyMDENuTVSppWXlIiIhA7cT05GRrYa9jQqONkrYW6vgYK2Cg40SDjYqONioYG+lhEIhyB0yPYnabQpeP7hjmnte3wOc+lP3GEBOBuBcCwjoC3R+Q/e6LEm3gavbgIijQNxlIP2+rt3BA6jTFmg5FmhQ9BErAMAf3YDYUMA9AHjjTOn3yUgCfmwIaHJ0jzQM/Kl83ycRmT2DklAvLy9ERkbiypUr+SOcAPDff/8BALp06VLknPT0dAC6hUpEVVl6thoRiemISMhARGI6biekIyIhHRGJ6Uh4mFOua9lZ6ZJSDycbtPF1RXv/6mjv5wYvF9sKip6MQmlV8FrUVPz9dn0MnJwnbUu6qWsLXQOMWVf6+ckRwC8tiz+Wt9jq0kag+Qjg+d+LrvpvPR7YMR1ICAfunNE9mlCS0LW6BBQAWo0tPS4iskgGJaEdO3ZEREQE5s+fj7Fjx8Le3h63bt3C5s2bIQgCnnnmmSLnhIeHA9AlsERVSUpGLnZfjsWusFiERacgLjW72H4eSEag8KB8F1cD93NdceVhNVyJScXyk1EAAB83e7T3d0N7Pze093eDb3V7CAJHTSuNuEsFr50q+D3txO8FCahTTd2jALXbAOos4Pp/wMn5wLoJQG4ps0xaDaC0Buo9DdR7CvBoqHusIDMZSLwBnF4AxF/RJbTV/ICnPpae33w4sPszXYWA4OWlJ6HBy3WfazQFarM2NBEVZVASOmXKFKxevRqhoaFo2rQpWrdujcOHDyMrKwv29vYYPXp0kXMOHz4MAAgICDDk1kQmkZKZiz2X47A9NBpHbyQgV1P2Cugxqn14R/Vvue81Rz0Ec9RBkraopAxEJWVg/bm7AAAPJ5v8pLRzvepoUMOp3PchI9GogRO/FXzt163kvoZ6GA/s/1r32sUHmLJXWhLKr4susVw+BNCqS76OkxfwzsXiE+a6PXXT5ptfB4JXAMd/Azq9Dtg+VlDf1gVo8jwQuhoI2wj0/R6wsit6rZgQIPai7nXLMeX+donIMhiUhPbq1Qtvv/025s6di4iICERGRuY/Jzpr1iy4u7tL+mdlZeWPknbv3t2QWxNVmNSsXOy9HIftoTE4fD1er8TzcSvUT2OPpg0EAfB0soGLnRXEnEzMyfgfAGBE7hd4qLUuct590bXMa8enZWN7aAy2h+oqU7TxrYYJnf3Qr6kXrJTchdckctKB6AvAwe+Bu4+ei3Tx0RWtryghKwtGOJ/9pviapHV7AK0nlF5z09pB91ESQQD6fAOErAZy04FbB3VJ5+Naj9clodkpwJWtutHRwi6s0H1WWuum9omIimHwNh8///wznn76aaxbtw6xsbGoWbMmxo8fj169ehXpu2XLFjg7O8PFxQWDBg0y9NZERpOWlYu9Vx4lnuEJyNFo9T63lost/Nwd4OfuAP/qjz6728PbzR42KqWuU0468K0uCV3zxSvIUdghI0eNh9lqZORo8DBbjYdZalyOScXp20k4E5GEtKxSRrQeOReZjHORyajhbIOxHXwxqoMP3B1tnuhnQCU49L3uoyQOHsDIFYCqAn/utw7qPtu6Ag0HlNyv1djyFX7X5AIP7wM5D3VT9Xns3YD0eCA2rGgS6tcFqF5fN31/YXnRJFSdA1x89GxqQF/Aobr+8RCRRTHKXnMDBw7EwIElbHH3mOHDh2P48GJ+ayaSSUpmLhYcuYV/jt5Gek7ZC0usVQr0DPBA36ZeCKzlAt/q9rC1Upb7vtYqBaxV1nC1l46Idg/wwNQe9aDRirgWm4YzEUk4fTsJpyOSEJ9W/DOoABCXmo0f94Tj1/03MLBFTUzq7I9mdYrZl5yMx9VXl6B1fgtwrOCC63GXdZ9rNi99i1Cv5rrRR00pC+M0ucC5xbrRztjQ0vtmJBbf3mocsPcL4PZh4EFUQeF+ALi2A8hMKuhHRFQCvZPQuXPnYtCgQahbt25FxkNkEhk5aiw6FoE/D91EahkjjtZKBboHeGBg85p4urEnnGytSu1vDEqFgCa1nNGkljMmdPaDKIqISMzAmdtJOHU7CUeux+N+MUlpjkaLf8/fw7/n76GNbzVM7OyHvpyqN8zjW2oKgm7E07669FnJipaZrPtc1u5CSpVuodHDuOKPZyQBy14AYoL1u686q/j2lqOB/d8A2lwgeBXQ88OCYxceLUhyqgXUf1q/+xCRRdI7CX333Xcxbdo0NGrUCIMGDcLAgQPRpUsXrtSlKiUrV4OVp6Lw+8EbpZZR0iWe7hjQvCaeblwDziZIPEsjCAL83R3g7+6A4e28kavR4r9LsVhyPAJnIpKLPefxqfpxHX0xqYs/HGyMMvlhWfJ2TKoUDHy/3fW/ggS00UDd9H2NQN33qLLVJdkA8FMgkHq35K1IHT2BgGd19UaDVwA9PtCdmxoD3Nyv69NiJKAo/ywBEVkOvf9Hql+/Pm7cuIErV67g6tWrmDVrFtzc3NCvXz8MGjQIffv2hZMTV+pS5ZSr0WLDubv4Zd91RKcUP7qjEIAeAR4Y2LwWejepARc7eRPP0lgpFRjYvBYGNq+FsHspWHw8AltCopGjLvosa1xqNmbvDseKU1H4YlATPBvoxV8eqxo7V93oZl5h+ZJo1AWjpoVlpQJhj6o2NBsODP275OtkPSg7ptYTdEnog0hd4Xv/bkDIqoJ6qawNSkRl0DsJDQ8PR3h4OLZu3YqtW7fi2LFjSExMxPLly7FixQqoVCp069YNgwYN4rQ9VRparYitodH4eU84IhJLrp/Yv5kXpj0TgPqeVe8Xqaa1XTB7WAt81K8RVp+5g2UnIhGbWjTRjknJwtTl59GzoQe+ei4QvtVLWSVNlYtnE10SGntRl2iW9Fxo3MWSn/FMuqmbPgeApkNKvld8uG6hUlnq9wacawOp93Sjof7ddJ8BwKczUL1e2dcgIotWrgfFAgIC8N577+HgwYOIj4/HihUrMHLkSLi4uCA3Nxf79+/HtGnT0KBBAzRp0gQffvghjh49ml+2iciU9l6OQ7+5R/D26uASE9CnGnpg25td8fuYNlUyAX1cdUcbvP5UfRz58Cn8NroV2vkVvyvZwWvx6PPzYfyy7zqyck2wyw8Zrm5P3efMZCB8Z8n98p7HLM7jq99z0kvud/Yf/WJSKApqgF7eDFzfq1sxD3AUlIj08sSrFVxdXTFq1CisXLkS8fHx+Qlo/fr1IYoirl69itmzZ6NHjx7w9PTE+PHjsXbtWqSmphozfqIiUjJy8daqC5iy9CyuxaUV26eDvxvWT+2ERZPao2lt81pFnjdVv25qZ2x7syu61ncv0idbrcVPe8LRb+4RHA6PlyFKKpeWowHVo6Lw/32sK6tUWMRR3ar3krjVRf4zpSGrin/e89pO4PRf+sfVaqzumrkZwObXdG3WTkDgYP2vQUQWyyhLZpVKJXr27InZs2fj2rVr+c+Mdu/eHUqlMn/aftSoUfDw8EDv3r0xd+5c3Lx50xi3J8p37EYCnp1zGFtCoos93qKOC5a92B6rX+6Itn5uJo7O9JrW1n2/v45qBU+nonUsbyekY/w/p/H6yvOILeFZWaoEHD2BXp/oXj+IAv7sAZz+G7h3Dog8Duz9Elg2RLci3b7oLx0AdLU/G/TRvb6xF1g2GLi8RVd4//oeYPMbwOoxuu06S7pGYdV8C0Zp81bkBw4uvSA+EdEjFVK3JW/a/sCBA7h//36J0/YBAQEIDAzE1q1bKyIMsiBZuRrM2HoZYxacKvZ5yIY1nPDnuDbY9HoXdGvgYVELcwRBwKAWtbDvvR6Y3MUfimK+9e2hMXj6x4NYcOQW1OUo1E8m1PlNoMNU3eu0aGDHdODvXsCifsDRnwEbR2D4YsDKvuRrDPwJcPHWvb51EFg7DvirJ7AiCLiwDHCuBYxcWfo1Cms9Xvo1a4MSkZ4qvHhg4Wn7AwcOSKbtr1y5ggsXLlR0GGTGwu6lYNCvR/HPsdtFjjnbqjB7WAvseLubxa8Kd7K1wueDmmDrm13R2se1yPH0HA2+2X4FA389iquxfGymUuo3Exi9TrdPvF01XVklt7q65PSVI0DtNqWf71IHeOWwrsB+9fqA0gawcQFqNAN6/A+YegTwbFS+mBoNBKwdda/dAwCfDk/2vRGRxRFEGVcNhYeHY9u2bfD29sawYcPK7J+amgoXFxekpKTA2dnZBBFSZabRivjz8E38vCe82P3du9Z3x6xhzVHTxU6G6ArJSQe+raV7/XG07NOVWq2Idefu4LudV/EgI7fIcTsrJb4f2gzPt6wtQ3QW5MtHzyP3+B/w1EfyxvKkEm8Cv7bWve79FdD1nfKdf+C7gm1Rv0wxamhEJA998zVZK1cHBARg2rRpcoZAVdSdpAy8uyYYZyOL1kS0USnwUb9GGN/JD4ri5p4JCoWAEe188EwTL8zceRVrzt6RHM/M1eDt1cG4EPUAH/dvDGsVd1yqUOnxBVtzWtvrnsusKi4s031WqIAWo/Q7JzkCyHlUsSKdC+OILBW3T6EqRRRFrDt3F19tuVTsXu+BtZwxZ0RLNKhRtcstmYqbgzVmBjXH8HZ18MnGMFyNlVYTWHw8AmH3UjBvTGvUcLaVKUoLcHah7gMAfLsCk7bLG4++Mh8UrMhvNABwqqHfeZteByKPVlRURFRFGCUJTUhIwIoVK3DkyBHcunULaWlp0GhKrz8oCAJXx1O5pGXlYvq6EPx3qei+2AoBeK1nfbz1dAOO2j2BNr5u2PxGF3y19TJWnoqSHDsbmYwBvxzFvNGt0KFudZkipErjYTyQnQqkxQIHv3u0Q5MAdOWsFhGVj8FJ6KpVq/Dqq68iLU03gqLvI6aWvECEyi/6QSYmLz5TZKQOAHzc7PHT8BYWUXKpItmolPj2hWZo5e2KTzeFIfuxLUATHmZj9IJT+KhfI7zY1Z//fo2lKj4DuedzIGSltK3dFKBWS/2vUVVGeomoQhmUhO7fvx9jx47NTzx9fX3RvHlzuLq6QqHgaBQZx8W7KXhxyRncT8sucmxkO298OrAJHG34ZImxDGvrjcY1nTF1+TncTc7Mb9doRXyz/QouRD3AzKDm/JlbOqU1UM0faDMBaP+K3NEQURVk0Or4Pn36YO/evXB1dcWKFSvQr18/Y8ZWBFfHW569l+Pw5qoLyCy0vaSbgzW+H9IMfQK9ZIqsnCrZ6nh9PMjIwdurg3GomB2V6ns64o+xbVDf01GGyIiIqDLTN18zaLjyzJkzEAQBX331VYUnoGR5Fh+7jZeXnS2SgNb3dMTm17tUnQS0inK1t8aiie3w9tMNUHj2/cb9h3j+t6PYeTFGnuCIiKjKMygJ1Wp1z4x16dLFKMEQAbpp3y+3XMKXWy9DW2icvlPd6tgwtTO83cqxows9MYVCwLvPBOCfCe3gbCudfk/P0eDVFecx/yAXGBIRUfkZlITWq1cPAJCenm6UYIjSs9V4ZdlZLD4eUeRYUJs6WDK5PVzsrUwfmIV7qpEntr3ZDY1rFp1WmbnrKmbuuqr3okQiIiLAwCR05MiREEUR//33n7HiIQsWl5qFEX+dwN4r94sce++ZAMwKas7ySzLyqW6Pf1/tjKGt6xQ5Nv/gTXy2OQzawkPXREREJTDof/TXXnsNTZo0wZw5c3D27FljxUQW6GpsKl6Ydwxh96R7llsrFZgzoiXefLoBywJVAnbWSswe1hyfDmhc5Njyk1GYtjYYuRptMWcSERFJGZSEOjo6YseOHWjUqBG6d++OTz75BKGhocjKyjJWfGQBDoXHI2j+CUSnSP/euNpbYdmL7TG4Ffcvr0wEQcCUbnXxQ1BzFN4VdVNwNF5dfh5ZuaVvVkFERGRQiaY8ly5dQq9evZCQkKD/jQUBarW6XPdhiSbzs/NiDN5YdQGaQtO4vtXtsWhiO9T1MJMSQFWwRJM+dlyMwdurLyBXI/3z61S3Ov6e0Ja1RImILJBJSjQBwNy5c9GyZUskJCRAFMVyfZBl2xUWizeLSUDb+FbDxte6mE8Casb6N6uJBRPawdZK+lZy4lYixi44hQcZOTJFRkRElZ1BwxQ7duzAu+++CwBQKBTo1q0bWrRowR2TqEx7LsfhjZXnoS6UgA5sXhOzh7WArZVSpsiovHoEeGDZix0wedEZpGUXzG4E33mAEX+exLIX28PT2VbGCImIqDIyKAmdNWsWAKB27drYsWMHmjVrZpSgyLztuxKH11acK5KAju3ogxnPNYWi8IOGVOm183PDqpc7Yvw/p5GUXjD6eS0uDcP+PIHlL3ZgbVciIpIwaLgyNDQUgiBgxowZTEBJLweu3sery88XeYZwVHsmoFVd09ouWPtKJ3gVGvWMTMzAsD9O4Mb9NJkiIyKiysigJFSj0a2AbdmypTFiITN3KDweryw/h5xCJXxGtPXG/w1mAmoO6ns6Yt3UTvCrLh31jE3NwvA/T+JqbGoJZxIRkaUxKAlt0KABACA5OdkowZD5Ono9AS8vPYsctTQBDWpTB98NacYE1Ix4u9lj7dROaFjDSdKelJ6DcQtPIyoxQ6bIiIioMjEoCR01ahREUcSmTZuMFA6Zo+M3EvDikjPILpSADmlVGzOHNmcCaoY8nWyx5pWOaOntKmmPT8vG2IWncD+NtYSJiCydQUnom2++ifbt2+PPP//E1q1bjRUTmZETNxMxuZgEdHDLWpg1rAWUTEDNlqu9NVZM6YB2ftUk7VFJGRi/8DRSMnNlioyIiCoDg4rVR0VFISUlBS+//DLOnDmDESNGYMSIEQgICIC9fdkrYX18fMp1Pxarr1pO3UrExEVnkFlo95xBLWrh5+EtoFJaUBkvMy1Wr4+UzFyM/OskrsRInwdt61sNy17sADtrluMiIjIn+uZrBiWhCoUifz9vURTLtbc3d0wyb2cjkjD+n9PIyJEmoAOa1cTckS0tKwEFLDoJBYD7aVkY9scJRBZ6HvSphh74a3xbWFna3wciIjNmsh2THt/9iDsmEQCE3UvBxEVniiSg/Zp6YY4lJqAETydbLH+xAzydbCTtB67FY/q6EGi1fD8gIrI0BhWrX7RokbHiIDNx70EmJi8+g4fZ0lHuPk1q4JdRrTjiZcG83eyx7MUOGPbHcaRmFfz92BwcDVc7K3z5XGC5ZlOIiKhqMygJnTBhgrHiIDOQmpWLyYvO4H5atqS9d+Ma+G10ayaghIZeTlg0qT3GLjgleVZ4yYlIVHOwxju9A2SMjoiITIlZARlFjlqLV5efw7U46a44HfzdMG9MK1ir+FeNdNr4VsMf49rASikd9Zyz9zoWH7stU1RERGRqBmUGBw4ceOJzX3vtNUNuTZWIKIr46N+LOHYjUdJez8MBf41rCxsVVz+TVI8AD/w0vCUKz75/ufUyNl24J09QRERkUgYloYMHD8a5c+fKfd7LL7+MP//805BbUyUyd991bDh/V9Lm7miDxZPaw8XeSqaoqLIb1KIWZjzftEj79HUhOHD1vgwRERGRKRmUhKalpaF///64du2a3udMmTIFCxYsMOS2VImsO3sHc/Zel7TZWSnxz8S28HYru1YsWbZxHX3x3jPS50DVWhFTl5/DhShuB0xEZM4MSkLr16+P+Ph4PPPMM7hz506Z/SdOnJi/on7kyJGG3JoqgaPXE/DRvxclbQoB+HVUKzSv4ypPUFTlvNGrPiZ38Ze0Zau1eHnZOcSkZMoUFRERVTSDktA9e/agdu3auHv3Lp555hncv1/8FJooihg3bhyWLl0KURQxduxYLFu2zJBbk8yuxqbi1eXnoC5U3/HL5wLRu0kNmaKiqkgQBHw6oDGGtKotaY9Py8ZLS88iI6d8m1oQEVHVYFAS6uvri927d6N69eq4fv06+vbti9RU6dZ8Wq0WY8aMwYoVKwDoyjotWbIECgVXS1dVcalZmLzoDNIK1QJ9qZs/xnfykycoqtIUCgEzg5qja313SXvYvVQWsyciMlMGZ4KNGzfGjh074ODggJCQEAwcOBBZWVkAAI1Gg1GjRmH16tUAgMmTJ+Off/5hQeoq7GG2GpMWnUF0SpakvX8zL3zUr7FMUZE5sFIqMG90a/i7S7c03XExFr/sv17CWUREVFUZZTiyXbt22LRpE6ytrXHs2DEMHToUmZmZGD58ONatWwcAeOmll7BgwQImoFWYWqPF6yvO43KMdLS7tY8rfhreEgoF/2zJMC72Vvh7fFs42Ur30Ziz9zq2h8bIFBUREVUEo82J9+rVC6tWrYJCocCuXbvg7++PTZs2AQCmTp3KkkxVnCiK+GzzJRwKj5e0+1W3x4IJ7WBrxVqgZBz1PR0xb3RrFP6d5r11wQi7lyJPUEREZHRGfTBz8ODB+PvvvwEA9+/fhyiKeP311/H7778b8zYkg6UnIrHqdJSkrZq9FRZPag83B2uZoiJz1T3AA58NbCJpy8rV4qWlZ3E/NauEs4iIqCrRa+/4qKiosjs90qtXL7z11luYO3cugoKC8P7775d4vo+Pj97XJfmciUjC19suS9psVAosmNAOfoWe3yMylomd/RAel4ZVpwvKv8WkZOHlZeew+uWOHH0nIqriBFEUy1x2qlQa/81eEASo1eUrvZKamgoXFxekpKTA2dnZ6DFRUXGpWRj461HEp2VL2ueNbo0BzWvKFFUVlJMOfFtL9/rjaMCaybs+ctRajF14CqdvJ0naX2hVGz8Nb8FnzImIKiF98zW9puNFUayQD6rcctRavLr8XJEE9LWe9ZiAkklYqxT4Y2wb1KlmJ2nfeOEe/jh0S6aoiIjIGPSajs/b5Ygsy1dbL+F81ANJW7cG7nivT0N5AiKL5OZgjYUT2mHI78eQnqPJb//hv6uo7+mIZ7g5AhFRlaRXEjphwoSKjoMqmbVn7mDFKemzvHWq2eGXka2gZCkm83XgO+DQ90XbJ2wD/LuV3PfLil213tDLCXNHtsJLy84ibxJFFIF3Vl/Ahtc6o5GXzI/n3D4CLBlYtL3H/4CnPjJ9PEREVQC3LaIiQu8+wKebwyRttlYK/DmuDapxJbzp3T4CfOmi+zjwXel90+KA39oX9N/2LmAmj770blIDH/ZtJGlLz9FgypKzSErPKd/F1NnAt7V1P6MrW40YJRER6UuvkVCyHAkPszF12TnkqLWS9u+HNEdgLReZojIzcWGAdwfjXzc1GlgyCEi8ofu6w6tAv2JGNfX16omC19V8DYvNSF7pXhfhcWn49/y9/La7yZl4Z00wFk9sp/+GCRFHgJyHgNIGqPuU4YHVbi39ec3vZPg1iYjMHEdCKZ9ao8UbK88X2ZJzUhc/DG5VW6aozFD8NeNf88EdYFH/ggS081uGJaAAUKNJwUclWc0vCAK+faEZWvu4StoPh8dj3oEb+l/o2i7dZ/9ugI2j4YFZO0h/XkREVCa9ktClS5di6dKlSE1NLbtzOaWkpORfn+T1/c6rOHlLWgqnvb8bPu7PPeErteQIYHF/IPm27utu04E+X8saUkWytVLij3Ft4OFkI2n/eW84jt9I0O8i4Y+S0IC+Ro6OiIj0pVcSOnHiREyaNAl37941egB3797FxIkTMXnyZKNfm/S3JSQaC47elrR5Odti3ujWsFJywLzSSrwJLBoAPHi0iKznx8DTn8kbkwl4Otni11GtJFt7akXgrdUXEFfWjkqxYUDKowL4TEKJiGRTaZ4JZd1Q+VyJScWH60MlbdZKBeaPbV1ktIkqkYTrumdA02J0Xz/9BdBtmrwxFRYTCiwfAqTHA441gLH/Al5NdccKr67PSgVOzAOubAGSIwErO6BWK6D7+4DPY8/QPowHTv2Bjle34ZpdBB5qVDirbYg56qG49NAPb668gJUvdYCqpF+e8kZBazQFXL2LHr+yFQheBURfADISAKU1YF8dcK4F+HUFAvoBddoY72dERGShypWExsTEwNHRCM9PPSY6Otqo16PyScnIxSvLziEzVyNp/+r5QLTyqSZTVFSm+1eBpc8BD+N0Xz/7LdDpdXljKizyOLByJJCdArj6AOM3A251i++bchdY+nzBM60AkJsO3NgD3NwPBC0EAl/QjWKuGAak6d43rABUE4BnlOfQXRGKibkf4EREIGbtvoaP+pXwGElJU/FaDbB+MnB5k7Rdk6NbxPQgEog6AVzfA7xyqNw/DiIikipXEtqnT5+KioNkoNWKeGfNBUQlZUjaR7X3xqj2PjJFRWWKDdMlbBkJAASg/yyg/UtyRyUV/h+wdgKgzgQ8GgPjNgLOpeyytXaCbnV/12lA/d66UdCok8DB74DsVGDzm7pR0ZUjdNfs9ZluVFJhhcwr/0F57EfYCLmYZfUnnsr+CX8euoW2vm5FC9k/jAfundO9bthPeuzMwoIE1KcT0Ho8UM0fsLYHMpKAuEvAjb26eIiIyGB6J6GcLjc/fx6+hQPX4iVtLbxd8eVzgTJFRGWKCwNO/wVkJgEQgEFzgDYTZQ6qkNB1wKapgFYN1G4DjFkP2LuVfk7sRWDSDqBO24K22q2B6vWAlcOBnDTg76cBiMBL+yUjqnZ12uCO1hHeJz5DHSEBvRQX8J+2Hd5bG4ztb3WDt5t9wTWv/weIWsDBQxfb4y5tfHTftrri/MpCb4/1ngI6v6FLSImIyGDcttNCnYtMxuzd0lJB7o7W+GNsa9iolDJFRWW6uq3g9YAfK18CevpvYMf7AETAvwcwcqV+JZA6vipNQPMEPAu4+AApUbqR3wE/FTul793rJahPfQOVNhvtFFfxn7YdUrPUeH3leayb2qng73TeVHyDZwGhUE3RvEcbvDsUTUAfV1ZCTUREeuG2nRYoJSMXb626AI22YHRbEIBfRrVCTRc7GSOjsgkAHv25he8CWo0DVJVkF6tDs4AD3+heNxoIBP0DqPRc2NZ0aMnHagTqklAIQNMhxfexsoPSvR5w/zJ8hPv5zaF3U/DNtiv4enBT3S5JNw/oDjQsZlW8kxeQdBMI3wl0ew9wqK5f7ERE9ERYe8fCiKKIDzeE4t6DTEn7m70aoHM9d5miIr21exHweLR15fXdwIbJgEYtb0wAsOujggS05Vhg+FL9E1AAqF6/5GO2j3bqsq8O2JW8WE6wdQUAuFvnStqXnYzElpDosndJajFK9znpFvBLK2DT68DF9UDKvaJ9iYjIYExCLczyU1HYdSlW0tbe3w1v9SolCaDKw95dusr8ylZg06uAVlv6eRXt5O+6z55NgOd+BRTlfKTD2r7kY8KjtymrUvo81q+Rpz1sraRvbf/bEIoHwY/2iPfrWvwjAq3H6UZAFSrdiv7g5cCGF4GfmwBzWwL/fQIk3S56HhERPREmoRbkcnQqvt52WdJWzd4Kc0e2LLmmIlU+Tl7AhK26ZyUB4OJaYNvbgJyLBxs/p/t8/zKw60P54gBgb63C1883lbRl5GiQdWmH7ovCq+If9/TnwFsXdKvv/bsXJL7Jt4ETvwG/tdOtoiciIoMx87AQGTlqvLHqPHLU0hGz2cNa8DnQqsilDjBhM+D0qOzR+aXAThmTv6B/dM+BArrV+7s+ki8WAMPaemN42zr5XzcUouAl6p4VFQOeLf1kVx+g+3Rdov+/KGDybqDDq4DKFtDmAtvfA2JCKjJ8IiKLwCTUQny++RJuxadL2iZ38cfTjWuUcAZVGI+GxrmOW11g/BZduSEAOP0nsOcL41y7vJRWQNCiggLwJ38Hdn8qTyyPzHi+KRp5OQEAnlacBwBc0fpgw81yvO0prXS7NfX7Hhi64FGjCFzebORoiYgsD5NQC7Dxwl2sP3dX0tastgs+7GekZIjKp0bTsvvoyyMAGLepYMHOsTnAwZnGu355qKyB4cuABo82tTj+q3xJMQBbKyXmj20DRxsVeit1SehebWt8sTkMEQnpZZxdDP8eBa8zEo0UJRGR5WISauZuJ6Tj041hkjZHGxV+HdWK9UDNhVdT3Z7sNs66rw9+q0sA5aCyBkYsB+o9rfv62Bxg39fyxALA390BPw+ojRbCTQDAfk0rpOdo8PaaYORqCi3mCllTeqWBm/sLXrv6VkC0RESWhUmoGctWa/DGyvNIz5HuC/9/LzSFn7uDTFFRhajdGhizDrB69Oe6+1PgzILSz6koKhtdkfq8MkhHZgMHvpUnFgDPWIVAKYiIF50RLNYDAITceYC5e69LO258GfipMbBtmi4hvXNa9+zn9b26lfEbp+r6WTsCzYeb+LsgIjI/5do7nqqW73ZcxaVo6T7XI9p64/mWtWWKiCqUT0dg1CrdNpfqLGD7dEBlB7QaY/pYrGwLYrl9GDg0E1BYAT3eN30s4TsBAGet2kHMLvi9e97BG+jawB0d6z5WlD79PnB2oe6jODYuukVYLnWKP05ERHrjSKiZ2nM5DouPR0ja6ns64ovnmsgTEJlG3R666XClNQAR2PIGELZBnlis7IBRawDfrrqvD3wDHPnRtDE8tktS457DoVIUbNUpisC7a4KRkvGouP1rJ4HeXwEB/XQbAti5AYJSVyy/Tjugx/+AN88CDXqb9nsgIjJTgig+eXHB2bNnY/z48fD09DRmTCVKTU2Fi4sLUlJS4OzsbJJ7VkXRDzLR/5cjeJBRsHOMjUqBLW90RcNHq4XJxHLSgW9r6V5/HA1YV9LHIQ58Bxz6Xvf6yxR5YzGGG3uB5UN1uyR9cAvzT8Rh5q6rki4DmtXEb6NbQSi8l7whvny0y1OP/wFPyVuuiojI1PTN1wwaCf3ggw/g7e2NF154AVu3boVW7l1bCBqtiLdXX5AkoADwxaBAJqBUPnGXCz5ynmA1eWUQ/p/u86Ndkl7pXhed6kr3hN9+MQbrClWPKLecdOnPi4iIymTwdHxubi62bNmCwYMHo3bt2vjwww9x9erVsk+kCvHHoZs4E5EsaRvQrCZGtfeWKSKqsuZ3Kvi4d17uaJ6MZxPdaGS3aQAAhULATyNawNXeStLtyy2XcPtJyjbluXde+vMiIqIyGZSEXrx4Ee+88w7c3d0hiiLi4uIwe/ZsBAYGonPnzli4cCEePnxorFipDKF3H+DnPeGStjrV7PDd0GbGnWokqiraTtJNh/t1zW+q6WKH74c0l3TLyNHg7dUXiuwoRkREFcegZ0LzqNVqbNu2DYsWLcLOnTuhVqvzkx57e3sEBQVh0qRJ6N69u0H34TOhJcvM0WDAr0ckuyIpBGDd1M5o41tNxsgIQNV5JtSCfPRvKFadviNpe7VnPXzYt5FMERERmQeTPBOaR6VSYfDgwdi8eTPu3r2LH374AY0bN4YoikhPT8fSpUvx1FNPISAgAN999x2io6ONcVt6zHc7rxTZlvONp+ozASUqwWcDm6Cuh/SXgT8O3cTxmwkyRUREZFmMXqLJ09MT06dPR1hYGE6ePImXX34Zzs7OEEURN27cwKeffgpfX1/0798fGzZsQG5ubtkXpVIdvHYfS09EStpa1HHBm083kCkiosrP3lqFX0a2gpVSWrZp2poQJKfnyBgZEZFlqNA6oe3bt8cff/yBmJgYLF26FF5eXhBFERqNBv/99x+GDx+O2rVr43//+x9iY2MrMhSzlZSeg/fXh0rabK0U+GlES1gpWQaWqDRNa7vg/WcbStpiU7Pw0b8XYYQnlYiIqBQVnqVERkZi5syZ+PzzzxEXF5f/rKgoihBFEQkJCZg1axbq1auHn3/+uaLDMSuiKOLjfy8iPi1b0v7JgCao5+EoU1REVcuUrnXRtb67pG3XpVisPnOnhDOIiMgYKiQJzcrKwvLly/H000+jXr16mDFjBiIiIiCKIho0aICZM2ciOjoau3fvxogRI6BUKpGZmYnp06dj+fLlFRGSWVp/7i52XZKOID/V0ANjO/jIFBFR1aNQCPhxeAtUK1S2acbWy7gVz+oeREQVxahJ6IkTJ/Dyyy/Dy8sLEyZMwMGDB6HVamFra4tx48bh0KFDuHr1Kt5//314eXmhd+/eWLVqFa5cuYIWLVpAFEWOhurpTlIGvtoqLYrt5mCNmUHNWY6JqJxqONti5lBp2abMXA3eWxcCtYZlm4iIKoLBSWhMTAxmzpyJxo0bo2vXrli4cCFSU1MhiiJatmyJefPmISYmBkuWLEG3bt2KvUa9evUwc+ZMAEB4eHixfaiARivi3TXBeJitlrR/N6QZPJ1sZYqKqGrrE+iFMYVmES5EPcCfh2/JFBERkXlTGXJy//79sWfPHmi12vyH+F1cXDB69GhMmTIFrVq10vtadevWBQBkZGQYEpJF+OPQTZyNlO6KNLxtHTwb6CVTRETm4ZMBjXHsRgIiEgveh+bsDUfPhh4IrOUiY2RERObHoJHQXbt2QaPRQBRFdOvWDUuWLEFMTAzmzZtXrgQU0BW17969u8EF7c1d2L2UIrsi+bjZ4/NBgTJFRGQ+7K1V+HF4Cygee6IlVyNi2poQZKs18gVGRGSGDBoJ9fT0xIQJEzBlyhQ0aGBYTcpatWrh4MGDBl3D3GU+2lpQrS0oHaMQgJ9HtICjjUF/lET0SBtfN0ztUQ+/H7yZ33YtLg0/7QnHR/0ayxgZEZF5MShzuXv3LlQqJj+m8v3OK7hZaFek13rWRxtfN5kiIjJPb/dugP1X7+NqbFp+21+Hb6F34xpo58d/b0RExmDQdPzKlSuxdOlSpKam6n3Ow4cPsXTpUixdutSQW1ucQ+HxWFJoV6RmtV3wdm/uikRkbDYqJX4e0bLIbkrvrQ1BeqEFgURE9GQMSkInTpyISZMm4e7du3qfExcXh4kTJ2Ly5MmG3NqiPMjIwfvrQiRttlaKR/9JclckoorQuKYzpj0j3U0pKikD/7fjikwRERGZF9kyGG6Jp78vt1zC/UK7In3cvzHqe3JXJKKK9HL3umjjW03StvJUFA5cuy9TRERE5sPkSahGo1thymdJ9bMrLBabgqMlbd0DPDCuo69MERFZDqVCwI/DWsDOSilp/3B9KB5k5MgUFRGReTB5Enrt2jUAgJsbH+4vS+LDbHyy8aKkzclWhZlDm3FXJCIT8XN3wCcDpKvi76dl47PNl2SKiIjIPJRrOPLw4cPFtp85cwYJCQmlnpudnY2bN29i9uzZEAQBLVu2LM+tLY4oivhscxgS06WjLV8OCkRNFzuZoiKyTGM6+GD35TgcDo/Pb9saEo0+TWpgUItaMkZGRFR1lSsJ7dmzZ5EROFEUy7XISBRFCIKAV155pTy3tjhbQ2Ow42KspO2ZJjUwpHVtmSIislyCIOCHoc3R5+dDSM0qWB3/2eYwtPd3Qw1nbpdLRFRe5Z6OF0Ux/6O4trI+6tSpg3nz5mHw4MHG/D7Myv3ULHy2KUzSVs3eCt++wGl4Irl4udji68FNJW0PMnLx4YZQLrQkInoC5RoJPXDgQP5rURTRq1cvCIKAhQsXwt/fv8TzBEGAra0tatasCW9v7yeP1gKIooiP/r2IlMxcSfvXg5vCw8lGpqiICACea1ELuy/HYXtoTH7bwWvxWH3mDka195ExMiKiqqdcSWiPHj2KbW/fvj2aNGlilIAs3fpzd7HvqrT8y4DmNTGwOZ87I5KbIAj45vmmOH07CfGPlU37ettldK3vDm83exmjIyKqWgxaHX/79m3cunULAQEBxorHokU/yMSMrZclbe6O1vj6+aYlnEFEplbNwRo/DG0uacvI0eCD9aHQajktT0SkL4OSUF9fX/j6+rLmpxGIoogPN4QirdCWgN++0AxuDtYyRUVExXmqkSdGtpM+WnTiViKWn4os4QwiIiqMez5WEitPR+HIdWmZqyGtaqNPoJdMERFRaT4Z0Bi1XaXl0r7bcRVRiRkyRUREVLXoNYQ5Y8aM/Neff/55se1P4vFrWbI7SRn4v+3S/ahrONvgi0GBMkVERGVxsrXCzKHNMXbhqfy2zFwN3l8fglUvdYRCwUoWRESlEUQ9aosoFIr80kB5224Wbn8Sj19LH6mpqXBxcUFKSgqcnZ2f+L6ViVYrYtTfJ3HqdpKkffGkdujZ0FOmqMjoctKBbx8tLvs4GrB2kDceMpqPN17EylNRkrYvBzXBxC4lVwwhIjJn+uZrek/HF64NWrj9ST4IWHIiokgCOqq9NxNQoiri4/5Fp+Vn7rqGyMR0mSIiIqoa9EpCtVpt/kdJ7U/yYeluxT/EzF1XJW21Xe3wyQCWuyKqKhxtVPghSLpaPjNXg/fXcbU8EVFpuDBJJhqtiOnrQpCVK03GZw1rDkcbVhsgqkq61HfH2I7SYvWnI5Kw+HiEPAEREVUBTEJlsujYbZyPeiBpm9DJF53rucsTEBEZ5KN+jVGnmnRa/of/ruJ2AqfliYiKwyRUBjfjH2LWf9ckbb7V7fFhv0YyRUREhnIoZlo+K1eL99eFQMNpeSKiIgxKQmNjYzF58mRMnjwZ9+7dK7P/vXv3MHnyZLz44otISkoqs7850mhFfLA+FNnqgml4QQBmBbWAvTWn4Ymqss713DG+k6+k7WxkMhYduy1TRERElZdBSeiyZcuwePFiBAcHo3bt2mX2r127NoKDg7F48WIsX77ckFtXWYuO3ca5yGRJ28TOfmjv7yZTRERkTB/2bQRvN+m0/Kz/ruFW/EOZIiIiqpwMSkJ3794NQRAQFBSk9zkjRoyAKIrYuXOnIbeukm6VMA3//rMNZYqIiIzNwUaFWUEtJG3Zai3eXx/KaXkioscYlISGhYUBANq3b6/3OW3btgUAhIaGGnLrKqe4aXgA+GFoc07DE5mZjnWrY2JnP0nbuchk/HOU0/JERHkMSkITExMBAB4eHnqf4+7uLjnXUiw+HoGzxUzDd6hbXaaIiKgifdC3IXyr20vaZu++hpucliciAmBgEuro6AgASElJ0fuc1NRUAIC1tbUht65SbiekY9Z/0qL0Pm72+KAvp+GJzJW9tW5a/vGdjbPVWkznankiIgAGJqF16tQBAJw4cULvc44dOwYAei1kMgdarYgP1hctSv9DEKfhicxde3+3ItPyF6IecFqeiAgGJqE9e/aEKIr49ddf80c4S5OamorffvsNgiCgZ8+ehty6ylh8PAJnIopOw3fkNDyRRfjg2Ubw47Q8EVERBiWhr7zyCgRBQExMDAYMGIC4uLgS+8bGxmLAgAGIjo6GIAh45ZVXDLl1lRCRkI4fOA1PZNHsrJX4oZhp+Q+4Wp6ILJxB88GBgYF4++23MWfOHBw/fhz169fHiBEj0K1bN9SsWRMAEBMTg8OHD2Pt2rXIyMiAIAh4/fXX0bJlS2PEX2lpH62GLzwNP5Or4YksTt60/KJjEflt5x4VsZ/Sra58gRERycjgbGj27NlISUnBokWLkJ6ejkWLFmHRokVF+omi7jf+KVOmYM6cOYbettJbciICpyOku0JN6OSLTvU4DU9kid5/tiH2X72PyMSM/LZZ/13D041rwN/dQcbIiIjkYfDe8QqFAgsXLsSmTZvQqVMnALqE8/EPAOjSpQu2bNmCv/76C8Lj81JmKCIhHTN3Safhvd3s8EFf7g1PZKnsrVX4Yah0b/lsNfeWJyLLZbR54eeeew7PPfcckpKSEBwcjISEBAC6uqCtWrVCtWrVjHWrSk2rFfHBhqLT8D8MbQEHG07DE1myDo+K2C8+HpHfdjYyGUuOR2ByV3/5AiOiSi0zR4PXVpzD60/VR1s/89nm2+hZkZubG3r16mXsy0qkp6dDqVQWaVcqlbC1tZX0K4lCoYCdnd0T9c3IyMgf4S1sxZl7OH1bOg0/qm0tNPeyLXIPQRBgb1+wajYzMxNarTR5fZyDg8MT9c3KyoJGozFKX3t7+/yR7OzsbKjVaqP0tbOzg0KhG5jPyclBbm6uUfra2trm/10pT9/c3Fzk5OSU2NfGxgYqlUr/vo99nZ6eAZQQhrW1NaysrAAAarUa2dnZJV738b4ajQZZWVkl9rWyssqvzVuevlqtFpmZmUbpq1KpYGNjA0A3W5KRkWGUvuX5d18Z3iPe6O6D/VfvIyqp4Hv6YddVdPR1hK+bdBU93yMKmP17xKO+5fl3z/cI83yPKPzvPiMjA9//dx0HrsXj4LV4jOtQB2895Q87K2WlfY8o7WchIVYhKSkpIoASP/r37y/pb29vX2LfHj16SPq6u7uX2Ldt27aSvr6+vsX2U7l6ib7v/Sv6frgt/8PvzaWiYGVbbH9fX1/Jddu2bVtiDO7u7pK+PXr0KLGvvb29pG///v1L/bk9LigoqNS+Dx8+zO87YcKEUvvev38/v+9rr71Wat/bt2/n950+fXqpfcPCwvL7fvHFF6X2PX36dH7fH374odS+Bw4cyO/722+/ldp327Zt+X0XLVokOeblKIitvBT5H7uXzBLFyOOi+IWzKH7hLHauo5Acz/vwchTERYsW5V9327Ztpcbw22+/5fc9cOBAqX1/+OGH/L6nT58ute8XX3yR3zcsLKzUvtOnT8/ve/v27VL7vvbaa/l979+/X2rfCRMm5Pd9+PBhqX2DgoIkf4dL6yv3ewQAsUmTJuLxGwmS9wnfD7eJNUZ9JwKCpC/fI8zzPaLwx9q1a/P7rl27ttS+fI/QfZj7e8TjGnbpJ/p8sEXyfuE+UPdvoLK/R6SkpIilMepIaFxcHA4ePIiwsDAkJelGA93c3NC0aVP07NkTNWrUMObtKh2HZr0BlXQnKJvgdRBzS/6NkszPK22s8WVPm4KGW18Dtwq+PPaiY7HnfXmw5BENMi+d6lXH+E6+WHoiMr/N1qcZnFoPQNr5bTJGRkSVSVauBtkth0MQCpbwiOocPDi+WsaojEd49FuBQWJiYjBt2jT8+++/JU6nqFQqDB06FD/++GN++abySk1NhYuLC6Kjo+Hs7FzkuNzD6KIoYsvFOHz33w2kZakxtqMPPu5TT+8h98oyjM6pNsOm2oSHcRDSC2rmWllZQ/XoumqNBrm5xV9XdKgBKzdvTrWVs29VnWpLz1bj2TmHcTe54GdmZ6XAxlfawbuanaRvHr5HmMd7RGl9OR3P94jH/93/3/bL+PvIbcnx956ui8mdfYr0BSrPe0RycjJq1aqFlJSUYvO1PAYnoSEhIejduzeSkpJK/IHm30wQUL16dezbtw/NmjUr973yktCyvim5xaZkYe6+cHw6oAkXIxFRiY7fSMDoBackbR3rumHllI5QKMy7iggRle5cZBKC/jiBx1Orlt6u2PBqZygr+fuDvvmaQSWa0tPTMWDAACQmJkIURfTu3Rtr1qxBREQEsrKykJWVhYiICKxduxZ9+vSBKIpISEjAgAEDSv2tparzcrHFd0OaMwElolJ1ru+OsR19JG0nbyVhxanIEs4gIkuQlavB++tCJQmotUqB2cOaV/oEtDwMSkJ/++03REdHQ6FQ4O+//8bu3bsxbNgw+Pj4wNraGtbW1vDx8UFQUBB27dqFBQsWQBAE3Lt3D/PmzTPW90BEVGX9r19j1Ha1k7R9t/Mq7iSZ7y/qRFS6n/aE41aCdHp/2jMBqO/pJFNEFcOgJHTz5s0QBAETJ07Eiy++WGb/yZMnY9KkSRBFERs3bjTk1kREZsHRRoUfgqRF7DNyNPhgfSi0LGJPZHHORSZjwZFbkrYW3q6YYoa1hA1KQsPDwwEAI0eO1PucUaNGSc4lIrJ0Xeq7Y3QH6bT8iVuJWHk6SqaIiEgOWbkavL8+BNpC0/A/DmsOldLgTS4rHYO+o4cPHwLQlWHSV97OSXoXMiUisgAf9WtUdFp+xxVOyxNZkJ/3hONWvDQ/ere3+U3D5zEoCfXw8AAAXLlyRe9zrl7V7anu7u5uyK2JiMyKk60Vvh8qrRqSnqPBhxs4LU9kCc5HJePvwtPwdVzwUjfzm4bPY1AS2rFjR4iiiJ9++qnU+m551Go1fvrpJwiCgI4dOxpyayIis9OtgQdGtfeWtB2/mYgVnJYnMmu61fCFpuGVCswe1sIsp+HzGPSdjR8/HgAQHByMAQMGIDo6usS+0dHRGDRoEM6fPw8AmDhxoiG3JiIySx/3L2a1PKfliczanL3XcbPQNPw7zzRAgxrmOQ2fx+Bi9UOGDMGmTZsgCAKsrKzQp08fdOjQAZ6enhAEAXFxcTh16hT27NmDnJwciKKIIUOGYP369eW+V1UpVk9EZIgj1+MxbuFpSRuL2BOZpwtRyRg6/7hkFLR5HRf8+2rnKjsKqm++ZnASmp2djfHjx2PdunW6CwrFv0Hm3WbYsGFYunRp/nZb5cEklIgsxccbL2LlKek0/IznAzG+k588ARGR0WXlajDw16O4cf9hfpu1UoFtb3VFQBUeBTXJjkmAbs/bNWvWYOvWrejXrx/s7OwgiqLkw87ODv369cO2bduwZs2aJ0pAiYgsSfHT8lcRlchpeSJzMXffdUkCCgBv925QpRPQ8jB4JLQwjUaDW7duISkpCYCufFPdunWhVCoNvjZHQonIkhy7kYAxhfaWb+/vhtUvcVqeqKoLufMAL/x+TDIN36y2Cza+VnWn4fPom68ZfXNzpVKJBg0aGPuyREQWp8ujveWXnyyYlj99OwlLT0RgYhfzLdtCZO6ycjWYXmg1vJVSMPvV8IVZzndKRFQFfdSvMepUk07Lz9x1DREJ3PCDqKr6eU84rheehn+6ARp6WcY0fB4moURElZhDMXvLZ+Zyb3miqupcZBL+KlSUvmltZ7zSo55MEclHr+n4yZMnG/3GgiBg4cKFRr8uEZG56VzPHeM6+mLZycj8ttMRSVh8PAKTu3JanqiqyMzRYPq6UIiFitL/NLwlrCxoGj6PXguTFApFiaWXnoQoihAEARqNplzncWESEVmq9Gw1+s49jDtJmflttlYK7Hy7O/zdHWSMjIj09eWWS1h8PELS9r9+jTDVzEZBjbowycfHx6hJKBERlY+DjQo/DG2BUX+fzG/LytXi/XUhWPNKJyi5Wp6oUjtxM7FIAtraxxUvdasrT0CVgF5JaERERAWHQUREZelUrzomdPLFkhMF0/JnI5Ox6NhtTLHg/8iIKruH2Wq8vz5E0mZrpdsb3pJ/gbS8BxCIiKqwD/s1go+bvaRt1n/XcCv+YQlnEJHc/m/7FdxNzpS0ffBsI9T1cJQposqBSSgRURVib63CrEKr5bPVWkxfFwINV8sTVTqHwuOx6rR0C94O/m6Y2NlPnoAqESahRERVTIe61Yv8B3Y+6gH+PHxTnoCIqFgpmbn4cH2opM3eWolZQS246xmMmITu27cP48aNQ/369eHo6AiVSoXLly9L+hw+fBi///47li9fbqzbEhFZpA/6NoRfdem0/M97wnE5OlWmiIiosBlbLyM2NUvS9smAxvAp9G/XUhmchGZkZGDYsGHo06cPVq5ciVu3biEjIwPFVX5SKpV44403MGHCBFy/ft3QWxMRWSx7axV+HN4Cjw+m5GpETFsbjGx1+crfEZHx7bkchw3n70raujVwx+j2PjJFVPkYnIQOHz4c//77L0RRRLt27TB9+vQS+3bp0gVNmzYFAGzYsMHQWxMRWbQ2vm5Fdlm5GpuGOXv5Sz6RnJLTc/DRvxclbU42Kswc2pwlLx9jUBK6YcMG7NixAwDw119/4eTJk/jhhx9KPWfIkCEQRRGHDh0y5NZERATgnd4N0KjQftN/HrqJsxFJMkVERJ9vuYSEh9nStkFNUMvVTqaIKieDktAlS5YAAMaOHYspU6bodU6bNm0AAFeuXDHk1kREBMBGpcTPI1rCSlkwuqIVgffWhSA9Wy1jZESWaXtoDLaGREvanm7kiaA2dWSKqPIyKAk9e/YsBEHAiBEj9D6nZs2aAID4+HhDbk1ERI80rumMac80lLRFJmbgu538ZZ/IlOLTsvHpJuk0vIudFb4b0ozT8MUwKAlNTEwEANSqVUv/Gyp0t9RqtYbcmoiIHvNy97po41tN0rb8ZBQOhfMXfiJTEEURH2+8iOSMXEn7jOcD4elsK1NUlZtBSaiLiwsAIDo6uoyeBW7fvg0AcHd3N+TWRET0GKVCwI/DWsDOSilp/2B9CFIK/adIRMa39uwd7LkcJ2nr19QLz7XQf6DO0hiUhAYEBAAAQkJCyuhZYNOmTQCAVq1aGXJrIiIqxM/dAZ8MaCxpi0vNxudbwmSKiMgyRCam46ut0tro1R2s8c3gppyGL4VBSeiAAQMgiiJ+/fVXZGVlldn/yJEjWL16NQRBwKBBgwy5NRERFWNMBx90D/CQtG0Ojsa2UP1nrIhIf2qNFu+uCUZGjrQ+7w9BzVHd0UamqKoGg5LQ119/HW5uboiLi0NQUBCSkoovCaJWq/H3339j4MCB0Gq18Pb2xsSJEw25NRERFUMQBPwwtDmcbVWS9k83heF+atmDBURUPvMP3sT5qAeStlHtffB04xryBFSFGJSEOjs7Y82aNVCpVNi5cye8vb3Rv3///OMffPAB+vTpA09PT0ydOhVpaWmwsbHB2rVrYWVlZXDwRERUlJeLLb4e3FTS9iAjF//792Kxu9kR0ZMJufMAc/dJN4fwq26PTws9FkPFM3jHpKeffhr79++Hj48PMjMzsWvXrvznH3bu3Il9+/bhwYMHEEUR3t7eOHDgANq3b29w4EREVLLnWtTCgOY1JW37r97HmjN3ZIqIyLxk5Kjx7ppgqLUFv9gpFQJ+GtESDjaqUs6kPHonoaUtPurSpQuuX7+OpUuXIigoCL6+vrCzs4O1tTVq1qyJAQMG4M8//8T169fRoUMHowROREQlEwQB3zzfFB5O0mfSvt52GVGJGTJFRWQ+vt1xBbcS0iVtrz9VH619qpVwBhUmiHrOzSgUCjRv3hzjx4/H6NGj4eXlVdGxFZGamgoXFxekpKTA2dnZ5PcnIqpq9l+Nw+TFZyVtbX2rYfXLHaFSGjwZRmSRDly9j0mLz0jaWtRxwfpXO8OK/670ztfK9ZO6ePEi3n//ffj4+KBfv35YvXq1XqviiYhIHr0a1cCo9t6StrORyfj94E2ZIiKq2hIfZuP99aGSNjurvO1zmYCWh94/rZEjR8LW1haiKEKtVmP37t0YM2YMvLy8MGXKFBw6dKgi4yQioif0yYAm8HGzl7TN3Xcd5yKTZYqIqGoSRREf/XsRCQ+zJe2fDGiMuh6OMkVVdemdhK5cuRJxcXFYtGgRevXqBUEQIIoiUlNT89v8/f3x+eefIzw8vCJjJiKicnC0UWHOyJZQKgqKZmu0It5efQGpWdxNiUhf687exe5CuyL1auSJMR18ZIqoatP7mdDC7t27h+XLl2P58uW4dOlSwQUfrYxv164dJkyYgJEjR6JaNeM8pMtnQomIntwv+67jpz3SQYLBLWthzkjuYEdUlqjEDPSbexjpjxWld3Owxq53usHTiXvDP07ffO2Jk9DHBQcHY+nSpVi9ejViY2N1F36UjFpZWaF///4YP348Bg4cCJXqycsWMAklInpyGq2IUX+dxOkI6cYic0a0xOBWtWWKiqjyU2u0GPHXySKPsPw1rg36BJp+oXZlVyELk0rSsmVL/PTTT7h79y527tyJ0aNHw87ODqIoIicnB5s3b8bQoUNRs2ZNvPnmmzh9+rQxbktEROWgVAj4eWRLOBWzmxLLNhGVbP7Bm0US0BFtvZmAGsgoI6HFSU9Px4YNG7Bs2TIcOHAAWq1Wd8NHI6QBAQG4cuVKua7JkVAiIsNtDYnGm6suSNpa+7hi7SudWLaJqJDQuw8w5PfjkqL0Pm722PF2NziyKH2xTDoSWhwHBweMHz8ee/bsQVRUFGbOnIn69etDFEWIosjFS0REMhnUohaC2tSRtJ2PeoBf9t+QKSKiyulhthpvrbogSUAVAvDziJZMQI2gwn/lVavVOHv2LE6fPo2oqKj8kVAiIpLPl88Fwre6tGzTb/uv40yh50WJLJUoivh040VEFHpU5fWn6qONL3dFMoYKS0JPnjyJ119/HTVr1sQLL7yAf//9Fzk5ORBFEY6OjpgwYUJF3ZqIiMrgaKPC3JGtoHqsbJNWBN5ZHYyUTJZtIlp37i42BUdL2lp4u+KtpxvIFJH5MWoSeuvWLXz11VcICAhAly5d8McffyAxMRGiKEIQBDzzzDNYtmwZ4uLi8M8//xjz1kREVE4tvV3x7jMBkrZ7DzLx6aYwVNByAaIq4cb9NHyx+ZKkzclGhd9GteKuSEZk8AMNycnJWL16NZYvX46TJ0/mt+e9gQUGBmL8+PEYO3YsatasaejtiIjIiKb2qIfD4fE4dbtgGn5rSDR6BnhgaKHnRoksQVauBm+svIDMXI2k/fuhzeFdaOcxMswTJaG5ubnYunUrli1bhp07dyI3Vzd1k5d4enp6YtSoURg/fjxatWIRZCKiykqpEPDziJboN/eIZBr+881haONbDX7uDjJGR2R632y/jKuxaZK20R18MKA5B9KMrVxJ6NGjR7Fs2TKsW7cOKSkpAAoSTxsbGwwaNAgTJkxA3759oVQqjR8tEREZXS1XO3w3pBleW3E+vy09R4O31wRj/dROnH4ki7HzYgyWn4yStDWs4YTPBzaRKSLzpncSWrduXURGRgKA5Fmhzp07Y/z48RgxYgRcXFyMHyEREVW4/s1qYkRbb6w5eye/LeTOA/y0Jxwf9m0kY2REpnEnKQMfbAiVtNlaKfDb6FawteLAWkXQOwmNiIjIf+3v749x48Zh/PjxqFu3bkXERUREJvb5oCY4HZGE2wnp+W3zD95EO79q6NWohoyREVWsXI0Wb62+gLQstaR9xnNN0aCGk0xRmT+951hcXFwwZcoUHD58GDdv3sSXX37JBJSIyIw42Kjwy8hWsFJK6zm/uyYEd5O5rSeZrx93h+NC1ANJ2/Mta2FYWy7Oq0h6J6GxsbH466+/0LVr14qMh4iIZNSsjgs+6tdY0paSmYvXV15AjlorU1REFedQeDz+OHRT0uZb3R7fDG7KDXYqmN5JqI2NTUXGQURElcSkLn7o19RL0hZy5wG+3XFFpoiIKsb91CxMWxMsabNSCvhtVGs42VrJE5QF4ZJHIiKSEAQBM4Oaw6/Qtp6Lj0dge2iMTFERGZdGK+KdNcFITM+RtH/UrzGa1eFCa1NgEkpEREU421rh9zFtYKOS/jfx4YZQ3Ip/KFNURMYz/+ANHL+ZKGnr3dgTk7r4yROQBWISSkRExWpSyxkzng+UtD3MVuO1FeeRmaMp4Syiyu9MRBJ+3ntd0lbTxRazglrwOVATYhJKREQlGt7WG0NbS1cIX41Nw+ebw2SKiMgwcalZeG3FeWi0BTXPFQIwd2QrVHOwljEyy8MklIiISiQIAr4Z3BQNC9VKXHfuLtY+VtieqCrIUWvx2orziE/LlrS/0zsA7f3dZIrKcjEJJSKiUtlZK/H72NZwsJbuGvPZpjBciUmVKSqi8pux7RLORSZL2ro1cMfrT9WXKSLLxiSUiIjKVM/DEd8PbS5py340qpSWlStTVET6W3v2TpF94etUs8MvI1tBqeBzoHJgEkpERHoZ1KIWJnTylbTdTkjH/zZchCiKJZxFJL/Quw/w6Sbpc8w2KgX+HNeGz4HKiEkoERHp7eMBjdGiUA3F7RdjsOR4hDwBEZUh8WE2pi47V2THr++HNkNgLdYDlZNKn05RUVFld3oCPj4+FXJdIiKqGDYqJeaNaY0BvxxFSmbBNPz/7biC5t6uaO1TTcboiKTUGi3eXHUB0SlZkvaJnf3wQivuCy83QdRjDkWpVJbVpfw3FgSo1epynZOamgoXFxekpKTA2dnZ6DEREZF+9l+Nw+TFZyVtHk422PJGF9R0sZMpKiKpb3dcwV+Hb0na2vu7YcWUDrBScjK4ouibr+n1JyCKYoV8EBFR1dSrUQ282rOepC0+LRsvLz3HQvZUKWwNiS6SgNZwtsG80a2ZgFYSek3HL1q0qKLjICKiKua9ZwIQevcBjt0o2Prw4r0UTF8fgt9GteLOMySbq7Gp+GB9qKTNSilg/tg28HCykSkqKkyv6fjKgtPxRESVy4OMHAyedwwRiRmS9mnPBOCtpxvIFBVZspSMXDw37ygiC/2d/PaFZhjdgWtRTMGo0/FERETFcbW3xoIJ7eBkK51Y+2lPOHZejJEpKrJUWq2Id9ZcKJKAjmznzQS0EmISSkREBqnv6YhfR7VC4Xrf09aGIOxeijxBkUWas+86DlyLl7S1qOOCL58LlCkiKg2TUCIiMljPhp74ZEATSVtmrgYvLz2L+2lZJZxFZDy7wmLwy77rkrbqDtaYP7YNbK2MX+WHDKfXwiR9iKKI4OBghISEICEhAZmZmWWugP/888+NdXsiIpLZ5C5+CI9Nw5qzd/LbolOy8Mqyc1j1UkcmAlRhzkcl4+3VwZI2pULAvDGtUcuVJcMqK6MsTFqyZAm++uorREZGlus8jaZ8ZTy4MImIqHLLUWsxdsEpnI5IkrQPaV0bPw5rwRXzZHSRiel44ffjSErPkbR/NrAJXuzqL1NUls1kC5M++eQTTJ48GREREXrVBWWdUCIi82WtUmD+2NaoXWj06d/z94rUbCQyVHJ6DiYuOlMkAR3R1huTu/jJExTpzaAk9NSpU/juu+8AAM888wyCg4Nx/vx5ALodkTQaDeLj47Fz504899xzEEURXbt2RUxMDLRabWmXJiKiKqq6ow0WTmwLB2vp9Pv3u65i35U4maIic5OVq8FLS8/idkK6pL17gAe+eaEpR92rAIOS0Pnz5wMAfH19sX37djRv3hxWVlb5xwVBQPXq1fHss89i06ZNmDdvHo4ePYq+ffsiJyenpMsSEVEV18jLGXNGtsLjeYAoAm+tuoBrsWnyBUZmQasV8d7aEJyNTJa0N/JywrzRrbgjUhVh0J/S8ePHIQgC3nrrLahUZa9xevXVVzF06FCEhobi999/N+TWRERUyT3TpAbef7ahpC09R4MpS88g8WG2TFGROZi56yq2F6pDW9PFFosntYeTrVUJZ1FlY1ASGhOj+wsQGFhQf0uhKLhkbm5ukXPGjRsHURSxZs0aQ25NRERVwKs96mFwy1qStjtJmZi0+AweZqtlioqqsmUnI/FnoeeLHW1U+GdiO3i52MoUFT0Jg5LQvCTT09Mzv83R0TH/dXx8fJFz6tSpAwC4ceOGIbcmIqIqQBAEfD+0OVp4u0raQ++m4OWlZ5GVW74qKWTZ9l2JwxebwyRtKoWA+WNbo3FNVs2pagxKQj08PADoluLnqVGjBpRK3cPoV65cKXJO3uhpWhqfCSIisgS2Vkr8Pa5NkRXzx28m4u3VF6DWcKEqlS307gO8sfICtIWK63w7pBm6NfCQJygyiEFJaN40/NWrV/PbrK2t89uLm3JftmwZAKBWrVpFjhERkXnydLbFshfbo7qDtaT9v0tx+GRjGMv2UanuJGVg8uKzyCw0cv7W0w0wvK23TFGRoQxKQrt16wZRFHHgwAFJ+4gRIyCKIv755x988cUXuHTpEk6fPo3XXnsNa9euhSAI6Nevn0GBExFR1VLXwxFLJreHo410Ieuas3fw/c6rJZxFli4lIxeTFp9BQqHFbENa18a7vRvIFBUZg0E7Jl26dAnNmjWDo6Mj7t69m18VPyMjA02bNkVERESROl2iKMLNzQ3BwcH5z4fqizsmERFVfSdvJWL8P6eRo5ZOw/+vXyNM7VFPpqioMspWazDhn9M4eUu6A1fnetWxeFJ7WKtYiqkyMsmOSYGBgThw4AA2btwItbpglaO9vT0OHDiALl26FNk1qWnTpti3b1+5E1AiIjIPHetWx7zRraFUSAcpvt95FatPR8kUFVU22WoNXlt+vkgCGlDDEfPHtmECagaMsnd8aa5du4ZLly5BrVajQYMGaNWq1RNfiyOhRETmY/25u5i+LkTSphCAeaNbo1+zmjJFRZVBjlqL11eex57L0h22PJxssPG1zqhTzV6myEgf+uZrZVeYN1DDhg3RsGHDsjsSEZFFCWpTBw8ycvDN9oJKKloReHt1MJxsrdC1gbuM0ZFccjVavLmqaALqaKPCoontmICaEY5lExGRbKZ0q4s3nqovacvRaPHysrMIvvNAnqBINrkaLd5efQH/XZImoA7WSiyZ3A5Na7vIFBlVBCahREQkq/f6BGB0Bx9JW0aOBpMWncaN+6wpbSnUGi3eXROMHRdjJe321kosntwebXzdZIqMKope0/FLly7Nfz1+/Phi25/E49ciIiLLJAgCvn6+KVIyc7E9tGA/8OSMXIxdcBqrX+4IP3cHGSOkiqbRinhvXQi2hUr3g7ezUmLRxHZo58cE1BzptTBJoVBAEAQIgiBZBZ/X/kQ3LnQtfXBhEhGR+cpRa/HikjM4cj1B0u7hZINlL7ZHIy++75sjjVbE9HUh2HjhnqTd1kqBfya2Q+d6fDa4qjF6iaa8EksltT/JBxERUR5rlQJ/jmuDVj6ukvb4tGyM+PMknxE1QxqtiA/WhxZJQG1UCiycwATU3Ok1HX/79u1ytRMRET0Je2vdCugxC07hUnRqfntKZi7G/H0Sf09oy8TETGi1Ij76NxQbzt+VtFurFPh7fFt0qc8/Z3NX4XVCjYnT8UREliE1KxcvLj6DMxHJknZrlQK/j26N3k1qyBQZGYNWK+KTTRex6vQdSbu1UoG/xrdBz4aeMkVGxmCSHZOIiIgqgrOtFZZO7oAeAR6S9hy1Fq8sP4fNwfdKOJMqO61WxGebw4okoFZKAX+Ma80E1IIwCSUiokrJzlqJv8e3Rf9mXpJ2jVbEO2uCsexkpEyR0ZPKytXg9ZXnseKUdHtWK6WA+WPaoFcjjnBbEiahRERUaVmrFPh1VGsMb1tH0i6KwGebwvD7wRsyRUblFZ+WjZF/ncTOMGkdUJVCwG98xMIi6bUwqVevXka/sSAI2Ldvn9GvS0RE5kWpEPD9kOZwtLHCP8ekC2J/2HUNaVlqfPBswycuGUgV73pcGiYtPoO7yZmSdpVCwK+jWuHZQK8SziRzplcSevDgQQiCUGpZpcL/+PP66ttORERUEoVCwGcDG8PZToU5e69Ljs0/eBNpWbmY8VxTKBT8v6WyOX4jAa8sP4e0LGltcCdbFf4c2waduQreYumVhHbv3r3UpDE6OhrXr+veFARBgJ+fH2rU0A2rx8XFISIiAqIoQhAENGjQALVq1TJC6EREZEkEQcA7vQPgZGuFr7ddlhxbfjIKaVlqzBzaHLZWSpkipMLWnb2Dj/69CLVWOohVp5odFk1shwY1nGSKjCoDg0s07dy5E2PGjIFWq8Unn3yCSZMmwd1d+ltNQkICFi1ahG+//RaCIGDFihXo169fue/FEk1ERAQAa8/cwf/+DUWh3AbNartg/tjWqFPNXp7ACIBu1vOnPeH4dX/RZ3ZbeLtiwfi28HCykSEyMgV98zWDktDw8HC0adMGKpUKR48eRWBgYKn9L1++jC5dukCj0eDs2bMICAgo1/2YhBIRUZ4dF2Pw9uoLyNVI/xtztbfCLyNboXuh8k5kGtlqDT5YH4rNwdFFjvUN9MLPI1rCzpqj1ebMJHVCf/zxR6Snp+ODDz4oMwEFgCZNmuCDDz7Aw4cPMXv2bENuTUREFq5/s5r4e3xb2BWafn+QkYsJi07jt/3XoS08VEoVKjk9B2MXnCo2AX25e138PqY1E1DKZ1ASumfPHgiCUK7V80899RQAYO/evYbcmoiICD0bemLT613g7+4gaRdFYPbucLy87CxSMnNlis6yRCSkY8j840V2uVIIwNeDm+Lj/o25cIwkDEpCY2Jiyn1O3gKn2NjYMnoSERGVraGXEza/0QV9iqkzuffKfTz321FciUkt5kwyls3B9/Dcb0dxOyFd0u5grcTCie0wrqOvTJFRZWZQEurq6goAOHTokN7nHDx4EADg4uJiyK2JiIjyOdta4c9xbfBh30YoPNgWmZiBF34/ho0X7soTnBlLycjFm6su4O3VwUgtVILJy9kW66Z2xlPchpNKYFAS2q1bN4iiiO+//x7h4eFl9g8PD8fMmTMhCAK6du1qyK2JiIgkBEHAqz3rYdmLHeDmYC05lpWrxbtrQvDF5jDkqLUyRWhejt1IwLNzDmNrSNHnP5vUdMam17ugSS0uIqaSGZSETps2DQqFAikpKejYsSPmzJmDpKSkIv2Sk5Mxd+5cdO7cGQ8ePIAgCHjvvfcMuTUREVGxutR3x7Y3u6KFt2uRY0tORGLkXycQm5Jl+sDMRFauBjO2XsaYBacQm1r05/h8y1pYO7UTvFxsZYiOqhKD64T+/PPPeO+99/Kf9RQEAf7+/vD09IQgCIiLi8Pt27chimL+bkmzZ8/GtGnTyn0vlmgiIiJ9Zat1ydKKU1FFjrnYWeGDvg0xqp0PF8uUQ9i9FLy7JhjX7z8scszZVoVvXmiG51pwQxpLZ5I6oXk2btyIN998E9HRBUPyeUnp45evWbMmfv31VwwZMuSJ7sMklIiIymvd2Tv4dFMYsouZhm/h7Yr/G9wUTWtznUJpNFoRfx6+iZ/3hBepywoAXepXx+xhLVDTxU6G6KiyMWkSCgC5ubnYvHkz9u7di4sXL+ZPy1erVg3NmjVD7969MXjwYFhZWT3xPZiEEhHRk7gUnYKpy8/hTlJmkWMKARjfyQ/T+gTA2fbJ/48yV3eSMvDummCcjUwucsxGpcD/+jXChE5+HFGmfCZPQk2BSSgRET2plIxcfLo5rNiFNADg7miDzwY2xnMtauXP5lmyHLUWq89EYebOq0jP0RQ5HljLGXNGtOT+71QEk1AiIqJiHL2egM83h+FWoZqWeTrVrY6vBweivqdlJlcarYiNF+5h7r7wEkeOp/aoh3d6B8BaZdD6ZjJTTEKJiIhKkK3W4O/Dt/Dr/hvFPitqpRTwUre6eLNXA4vZZlKrFbHrUix+2hOOG8UsPAIAbzc7/DS8Jdr5uZk4OqpKZEtCNRoNkpOTkZmZibIu7ePjU65rMwklIiJjupOUgS+2XML+q/eLPV7b1Q7Tnw3AgGa1zHbUTxRFHLwWj9m7r+FSdMk7Sw1vWwefDwqEo43KhNFRVWTSJDQhIQG//vorNm3ahMuXL0OrLbsQsCAIUKvVZfZ7HJNQIiIyNlEUsftyHL7acgnRJdQP9XCywej2PhjT0QeeTuZT//LkrUTM/u9asYuO8rTxrYbpfRqiU73qJoyMqjKTJaHHjx/HkCFDEB8fX+bIp+TGggCNpuiDzqVhEkpERBUlI0eNX/bdwIIjt6DWFv//mZVSwIBmNTGxiz9aFlMMv6oIufMAs3dfw5HrCSX2aVLTGe8/2xA9G3pwoRaVi0mS0MTERDRq1AiJiYlwdHTElClT4Orqii+//BKCIGDBggVISkrC2bNnsWXLFmRlZaFLly548cUXAQATJkyokG+KiIjoSV2PS8Onm8Jw6nbRHQAf19LbFRM7+6F/s5pVYqo++kEmdlyMwfaLMbgQ9aDEfvU8HDDtmYbo19SLZZfoiZgkCf3qq6/w1VdfwcbGBmfPnkVgYCAuXbqEZs2aFRnpjImJwejRo3H48GFMnz4dM2fOLPf9mIQSEZEpiKKI/y7F4Z9jt3G6jGTUw8kGYzr4YHSHyjdVH5OSiR0XY7E9NBrnS0k8AaBONTu80zsAg1vWgkpZ+ZNqqrxMkoR27NgRZ86cwdSpUzFv3jwAKDEJBYDMzEy0aNECN2/exJ49e9CrV69y3Y9JKBERmdrl6FQsOR6BTcH3il1Jn0elENDS2xXt/N3Q3t8NbXyryVL8PjYlK3/E81wpz3rm8XSywZu96mNEO58qMaJLlZ9JklB3d3ckJydj/fr1eOGFFwAAly9fRtOmTSEIAnJycqBUSktbzJ8/H6+//jqCgoKwdu3act2PSSgREcklKT0Hq89EYfmJyBIXMD1OIQCNazqjvb8b2vu5oZ2/G9wdbYwak0YrIvpBJiIS03EtNg3/XYrFmYiyE08AqO5gjVd61MW4jn4WU4aKTEPffM2gOgupqbpSDr6+vvlttrYFUxFpaWlwdXWVnNO2bVsAwKlTpwy5NRERkUm5OVjjtZ718XK3uthzOQ6LjkeUOlWvFYFL0am4FJ2KRcciAAB1PRzQwd8NzWq7wtlOBQdrFRxsVLC3VsLRRgV7GyUcrHVf5y0G0mpFxKZmISIhHbcT03WfEzIQkZiOqMQM5GjKrkiTx9lWhWcDvTCgeU10qe8OK067k4wMSkIdHR2RkpIiKbXk5lZQwDYiIgItW7aUnJOVpfvt8f794muyERERVWYqpQL9mtVEv2Y19Z6qz3MrPh234tOxCndK7ScIgL2VEvY2KqRm5up17ZI4PZ541nPnlDtVGgYlofXr18e5c+cQFRWF9u3bAwBcXV3h5eWFuLg4HDhwoEgSevToUQCAg4ODIbcmIiKSXZNazpgZ1ByfDWqCsxFJOH07CWcikhByJ6VcI5SFiSKQnqMpds92fTjZqPBMYA0MfDTiaaPidDtVPgYloR06dMC5c+dw5swZBAUF5bf37dsXixcvxg8//ICBAweiQYMGAICTJ09i1qxZEAQB7dq1MyxyIiKiSsLRRoWeDT3Rs6EnACArV4OQOw9w+nYSTkck4VxkMjKeMKHUh5ONCn7uDmjk5YRnA73QLYCJJ1V+Bi1M2rZtG5577jnUq1cP169fz28PCwtD69atodFooFQq0aJFC6Snp+P69evQaDQQBAHbt29H3759y3U/LkwiIqKqSK3R4nJMav5IaWxKlm6kM1ut+8jRQFNCgfw89tZK+FZ3gL+7PfzdHeBX3UH32d0B1R2sWVCeKg2TrI7Pzc3FSy+9BI1GgxkzZsDf3z//2MKFC/Hqq68WuzXnV199hc8++6zc92MSSkRE5kgURWSrtcjIS0xz1EjP1r22USng5+4ATycbJppUJZh07/iSXLt2DYsXL8alS5egVqvRoEEDjBs3Ln+FfHkxCSUiIiKq3CpFEmpsTEKJiIiIKjd98zXZ6jScO3dOrlsTERERkcxMnoQeP34c/fr1Q4cOHUx9ayIiIiKqJAwq0VQe+/btwzfffIPDhw+b6pZEREREVEmVOwkVRREbN27E3r17cefOHVhZWcHPzw9BQUHo3Llzkf4HDx7Exx9/nL9NZ94jqH369DEwdCIiIiKqqsqVhEZGRuL555/HxYsXixybO3cuhg0bhhUrVkCpVCIxMRFTpkzBli1bAOiST0EQ8Pzzz+OTTz554hXyRERERFT16Z2E5uTkYODAgbh06VKJfdatWwcfHx+8+eab6NGjByIjIyGKIpRKJYYPH46PP/4YgYGBRgmciIiIiKouvZPQFStW4NKlSxAEAb6+vvj000/RrFkzWFtb48qVK5g1axYuXLiA+fPn48SJE4iIiAAADB06FN9++23+1p1ERERERHrXCR00aBC2b98Ob29vXLp0CY6OjpLjWq0W3bt3x/HjxwEASqUSCxcuxPjx440WLOuEEhEREVVuRq8TGhISAkEQ8P777xdJQAFAoVBgxowZAABBEDBu3DijJqBEREREZD70TkITExMBAE2bNi2xT/PmzfNfBwUFGRAWEREREZkzvZPQzMxMAICnp2eJfdzd3fNf16lTx4CwiIiIiMicVdiOSSqVyergExEREVEVI9ve8URERERkuco9XPn777+XOiVfnn6ff/55eW9PRERERGZA7xJNCoUCgiAY9eYajaZc/VmiiYiIiKhy0zdfK9dIqJ75ql6MndASERERUdWhdxJ64MCBioyDiIiIiCyI3klojx49KjIOIiIiIrIgXB1PRERERCbHJJSIiIiITI5JKBERERGZHJNQIiIiIjI5JqFEREREZHJMQomIiIjI5JiEEhEREZHJMQklIiIiIpNjEkpEREREJscklIiIiIhMjkkoEREREZkck1AiIiIiMjkmoURERERkckxCiYiIiMjkmIQSERERkckxCSUiIiIik2MSSkREREQmxySUiIiIiEyOSSgRERERmRyTUCIiIiIyOSahRERERGRyTEKJiIiIyOSYhBIRERGRyTEJJSIiIiKTYxJKRERERCbHJJSIiIiITI5JKBERERGZHJNQIiIiIjI5JqFEREREZHJMQomIiIjI5JiEEhEREZHJMQklIiIiIpNjEkpEREREJscklIiIiIhMjkkoEREREZkck1AiIiIiMjkmoURERERkckxCiYiIiMjkmIQSERERkckxCSUiIiIik1PJHUB5iKIIAEhNTZU5EiIiIiIqTl6elpe3laRKJaFpaWkAAG9vb5kjISIiIqLSpKWlwcXFpcTjglhWmlqJaLVaREdHw8nJCYIgyB0OERERERUiiiLS0tJQq1YtKBQlP/lZpZJQIiIiIjIPXJhERERERCbHJJSIiIiITI5JKBERERGZHJNQIiIiIjI5JqFERAZYvHgxBEGAIAiIiIiQOxwioiqDSSgRWaSIiIj85NGQDyIiejJMQomIiIjI5FgnlIgsUm5uLq5du1bi8WbNmgEA2rZti0WLFpXYr2nTpkaPjYjIElSpbTuJiIzFyspKrwTSwcGBiSYRUQXgdDwRERERmRyTUCIiA5S1Or5nz54QBAE9e/YEANy4cQNTp05F3bp1YWdnBz8/P7z44ouIjIyUnBcWFoZJkyahbt26sLW1hbe3N1599VXcv39fr7g2bdqEYcOGwcfHB7a2tnB1dUXbtm3x1VdfITk52dBvm4jIYJyOJyIykb1792LIkCFIS0vLb4uMjMQ///yDbdu24dChQ2jUqBFWrVqFiRMnIicnJ7/f3bt38ccff2Dnzp04fvw4atWqVew9kpOTERQUhP3790vas7Ozce7cOZw7dw6///47Nm/ejI4dO1bMN0pEpAeOhBIRmUB0dDSGDx8OV1dX/Prrrzh16hSOHDmCd955B4Ig4P79+5gyZQrOnDmD8ePHo169eliwYAFOnz6NAwcOYNy4cQB0Seu0adOKvUd2djZ69+6N/fv3Q6lUYty4cVi1ahVOnjyJI0eO4P/+7/9QvXp13L9/H/379y8y+kpEZEocCSUiMoHr16+jQYMGOHbsGDw8PPLbu3btCpVKhdmzZ+PYsWMYMGAA2rdvjz179sDe3j6/X8+ePZGVlYV169Zhw4YNiI+Pl1wHAGbMmIHz58/D1dUVe/fuRZs2bSTHu3btijFjxqBTp06IiYnBxx9/jBUrVlTsN05EVAKOhBIRmcgvv/xSJHEEgNdeey3/dUJCAhYsWCBJQPO8+uqrAAC1Wo0TJ05Ijj18+BDz5s0DAHz99ddFEtA8vr6++OyzzwAA69atQ3p6+pN9M0REBmISSkRkAq6urnj22WeLPebv7w8nJycAQPPmzdG4ceNi+7Vo0SL/9a1btyTHDh06hJSUFABAUFBQqbF0794dgK5W6rlz5/T7BoiIjIzT8UREJtCgQYNSt/l0dXVFWloaAgICSu2T5/HFTQBw9uzZ/Nc1a9bUO67Y2Fi9+xIRGRNHQomITKC46fXHKRSKMvvl9QEAjUYjOaZv6abCMjIynug8IiJDcSSUiMgMPJ6Unj9/HlZWVnqdV6dOnYoKiYioVExCiYjMQPXq1fNfe3h4MLkkokqP0/FERGagVatW+a+PHTsmYyRERPphEkpEZAZ69+6d/zzpL7/8AlEUZY6IiKh0TEKJiMyAq6sr3njjDQDA8ePH8e6770Kr1ZbYPy4uDgsWLDBVeERERfCZUCIiMzFjxgwcOnQIp06dwty5c3Hw4EG89NJLaNmyJRwcHJCcnIxLly5h79692LlzJ5o1a4YpU6bIHTYRWSgmoUREZsLGxgZ79uzBxIkT8e+//yIkJCR/dLQ4zs7OJoyOiEiKSSgRkRlxcnLChg0bcPToUSxZsgRHjhxBdHQ0MjMz4ezsjHr16qF9+/YYMGAA+vTpI3e4RGTBBJFPrxMRERGRiXFhEhERERGZHJNQIiIiIjI5JqFEREREZHJMQomIiIjI5JiEEhEREZHJMQklIiIiIpNjEkpEREREJscklIiIiIhMjkkoEREREZkck1AiIiIiMjkmoURERERkckxCiYiIiMjkmIQSERERkckxCSUiIiIik/t/ApICpcMlaIsAAAAASUVORK5CYII="
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "id": "bfb66e0949b34052",
   "metadata": {},
   "source": [
    "## Companion Mass function\n",
    "\n",
    "如果我们在观测一颗恒星，我们发现这颗恒星的视向速度有变化，趋势如上:\n",
    "- 视向速度随时间的变化为正弦函数\n",
    "- 半振幅为$K$ [$\\mathrm{km/s}$]\n",
    "- 周期为$P$ [$\\mathrm{day}$]\n",
    "\n",
    "假设:\n",
    "- 可见星质量为 $m$\n",
    "- 伴星质量为 $M$\n",
    "- 以速度 $v$ 绕质心做圆周运动, 半径为 $R$ \n",
    "- 双星轨道平面倾角（平面法向量和视线方向的夹角）为$i$\n",
    "\n",
    "回答以下问题：\n",
    "1. 若已知这颗恒星的视向速度变化是由于其伴星引起的，求解伴星的质量的表达式？\n",
    "2. 如果我们发现这颗恒星的视向速度半振幅 $K=50 \\,\\mathrm{km/s}$ ，$P=1 \\,\\mathrm{day}$，倾角$i=0^\\circ$，求伴星质量？\n",
    "3. 如果视向速度半振幅 $K=52.8 \\,\\mathrm{km/s}$ ，$P=78.9 \\,\\mathrm{day}$，倾角$i=16.5^\\circ$ 呢？"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8d7072690ebacd2",
   "metadata": {},
   "source": [
    "## [Q0]\n",
    "\n",
    "假设可见星质量为 $m$ ，以速度 $v$ 绕质心做圆周运动, 半径为 $R$ ，伴星质量为 $M$\n",
    "- 万有引力定律：\n",
    "$\\dfrac{GmM}{R^2} = \\dfrac{mv^2}{R}$ \n",
    "- 周期、速度和半径的关系\n",
    "$P=\\dfrac{2\\pi R}{v}$\n",
    "\n",
    "伴星质量 $M = \\dfrac{v^2 R}{G} = \\dfrac{Pv^3}{2 \\pi G}$\n",
    "\n",
    "其中，$v$ 与视向速度曲线的半振幅 $K$ 之间存在关系：$v=K/\\sin{i}$, $i$为公转平面与视线方向的夹角（0度为face-on,90度为edge-on）\n",
    "\n",
    "因此， $M=\\dfrac{P K^3}{2 \\pi G {\\sin{i}}^3}$"
   ]
  },
  {
   "cell_type": "code",
   "id": "e3b03b4f6f78383f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-11-28T04:44:33.423957Z",
     "start_time": "2024-11-28T04:44:33.329018Z"
    }
   },
   "source": [
    "import astropy.units as u\n",
    "import astropy.constants as c\n",
    "\n",
    "const = (1/(2*np.pi) * u.day * (u.km /u.s)**3 / c.G).to(u.Msun)\n",
    "const"
   ],
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Quantity 1.03614907e-07 solMass>"
      ],
      "text/latex": "$1.0361491 \\times 10^{-7} \\; \\mathrm{M_{\\odot}}$"
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "id": "ec4481c722774956",
   "metadata": {},
   "source": [
    "因此，M的近似公式为\n",
    "\n",
    "$M =  \\dfrac{P[\\mathrm{day}] K[\\mathrm{km/s}]^3}{{\\sin{i}}^3} \\times 1.036 \\times 10^{-7} M_\\odot$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8eaf939dee4fa736",
   "metadata": {},
   "source": [
    "## [Q1]\n",
    "\n",
    "如果我们发现这颗恒星的视向速度半振幅 $K=50 \\,\\mathrm{km/s}$ ，$P=1 \\,\\mathrm{day}$，倾角$i=0^\\circ$，求伴星质量？\n",
    " \n",
    "$M = 1.036 \\times 10^{-7} \\times \\dfrac{1 \\times 50^3}{1^3} M_\\odot$"
   ]
  },
  {
   "cell_type": "code",
   "id": "b1fff323f26ee448",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-11-28T04:44:33.430487Z",
     "start_time": "2024-11-28T04:44:33.427672Z"
    }
   },
   "source": [
    "const *50**3"
   ],
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Quantity 0.01295186 solMass>"
      ],
      "text/latex": "$0.012951863 \\; \\mathrm{M_{\\odot}}$"
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "id": "720eecb679f937f",
   "metadata": {},
   "source": [
    "## [Q2]\n",
    "\n",
    "如果视向速度半振幅 $K=52.8 \\,\\mathrm{km/s}$ ，$P=78.9 \\,\\mathrm{day}$，倾角$i=16.5^\\circ$ 呢？\n",
    " \n",
    "$M = 1.036 \\times 10^{-7} \\times \\dfrac{78.9 \\times 52.8^3}{{\\sin{16.5^\\circ}}^3} M_\\odot$"
   ]
  },
  {
   "cell_type": "code",
   "id": "86930b18678c7023",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-11-28T04:44:33.441509Z",
     "start_time": "2024-11-28T04:44:33.438815Z"
    }
   },
   "source": [
    "const *78.9*52.8**3 / np.sin(np.deg2rad(16.5))**3"
   ],
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Quantity 52.52615332 solMass>"
      ],
      "text/latex": "$52.526153 \\; \\mathrm{M_{\\odot}}$"
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "execution_count": 5
  },
  {
   "cell_type": "code",
   "id": "f6201bb1872f6a7e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-11-28T04:44:33.457561Z",
     "start_time": "2024-11-28T04:44:33.455981Z"
    }
   },
   "source": [],
   "outputs": [],
   "execution_count": null
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
